/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAMOMDIS_M
  \brief Calculate residue momentum distribution with angular momentum component
  m specified. This is a class to assist class TAMOMDIS.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/09
  \date Last modified: 2020/10/06 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#ifndef _TAMOMDIS_M_h_
#define _TAMOMDIS_M_h_

#include <complex>
#include "TAMOMDIS.h"

typedef std::complex<double> cdouble;

class TAMOMDIS;

class TAMOMDIS_M{
public:
  TAMOMDIS_M(TAMOMDIS *mom);
  virtual ~TAMOMDIS_M();

  /// returning c.s. for m, mom is assigned to momStr
  double ParallelStr(int m, double *momStr);
  /// the same as ParallelStr, but for diffraction dissociation
  double ParallelDiff(int m, double *momDiff);
  double CSStr(int m) const; ///< total cross section for the stripping channel
  double CSDiff(int m) const; ///< total cross section for the diffraction dissociation channel
  double CSDiffEx(int m, int mp) const; ///< the left-out exchange term - 2023.1.2

  // accesses to the S-Matrices //
  TASMatrix *Sc() const{ return fSc; }
  TASMatrix *Sv() const{ return fSv; }
  cdouble Sc(double b) const;
  double Scr(double b) const;
  double Sci(double b) const;
  double Sc2(double b) const;
  cdouble Sv(double b) const;
  double Svr(double b) const;
  double Svi(double b) const;
  double Sv2(double b) const;

  const double *GetMomArr() const{ return fMomArr; }
  double kz(int i) const; ///< \retval fMomArr[i] in MeV/c

  static const int kNmom = 301; ///< number of points in momentum distribution

protected:
  double ParallelStr(int m, double kz) const; ///<\retval dsig_str/dkz
  double ParallelDiff(int m, double kz) const; ///<\retval dsig_diff/dkz
  double Rl(double r) const; ///< the radial wavefunction
  /// the sum_p{|B_lmp|^2} without z-integral. Ref[PRC.70.034609-12, the appendix]
  double Blm2(int m, double rho, double b) const;

  int l; // angular momemtum quantum number
  TAMOMDIS *fMOM;
  double *fMomArr; ///< momentum array for the momentum distribution
  /// parameters all from fMOM->fSc
  int fNG; ///< number of Guassians in Gaussian fit of the SMatrix
  double *fAlphajR, *fAlphajI, *fBetaM2, fRL; ///< fBetaM2 = 1/betaj^2, (beta^(minus 2))
  TASMatrix *fSc, *fSv;

  // set rho and b mesh //
  static const int NR = 301, NR2 = 301, NPHI = 101, NPHI2 = 101, NB = 301;
  static const double RMAX, BMAX, ZERO, MOMMAX;
  double bv[NB]{}, db; ///< bv mesh
  double r[NR]{}, r2[NR2]{}, dr, dr2; ///< rho mesh, 2: (-\infty,\infty)
  double phi[NPHI]{}, phi2[NPHI2]{}, dphi, dphi2; // 2: 2pi
  double fBLMSQ[NB][NR]{}; ///< integrate |sc|^2 over phi
  cdouble fSPINT[NB][NR]{}; ///< s'int, see Eq.5.22 of the DOC
  double fSINT2[NB][NR]{}; ///< sint, integral of |Sc|^2 over phi
  cdouble fSSINT[NB][NR]{}; ///< sint, integral of Sc*Sv over phi
  double fSSINT2[NB][NR]{}; ///< sint, integral of |Sc*Sv|^2 over phi
};

#endif
